home *** CD-ROM | disk | FTP | other *** search
/ Software Vault: The Diamond Collection / The Diamond Collection (Software Vault)(Digital Impact).ISO / cdr09 / vstsrc.zip / COMPUTE.C < prev    next >
C/C++ Source or Header  |  1995-01-23  |  29KB  |  871 lines
  1. /*
  2.  * %W% %E% %U%  [EXTREL_1.2]
  3.  *
  4.  * VersaTrack orbit calculations are based on those that in Dr. Manfred
  5.  * Bester's sattrack program (the Unix(tm) versions 1 & 2).
  6.  *
  7.  * The data from which the maps where generated come from "xsat", an
  8.  * X-Windows program by David A. Curry (N9MSW).
  9.  *
  10.  * Site coordinates come from various sources, including a couple of
  11.  * World Almanacs, and also from both of the programs mentioned above.
  12.  *
  13.  * The following are authors' applicable copyright notices:
  14.  *
  15.  *                                                                               
  16.  * Copyright (c) 1992, 1993, 1994 Manfred Bester. All Rights Reserved.        
  17.  *                                                                           
  18.  * Permission to use, copy, modify, and distribute this software and its      
  19.  * documentation for educational, research and non-profit purposes, without   
  20.  * fee, and without a written agreement is hereby granted, provided that the  
  21.  * above copyright notice and the following three paragraphs appear in all    
  22.  * copies.                                                                    
  23.  *                                                                              
  24.  * Permission to incorporate this software into commercial products may be    
  25.  * obtained from the author, Dr. Manfred Bester, 1636 M. L. King Jr. Way,     
  26.  * Berkeley, CA 94709, USA.                                                   
  27.  *                                                                             
  28.  * IN NO EVENT SHALL THE AUTHOR BE LIABLE TO ANY PARTY FOR DIRECT, INDIRECT,  
  29.  * SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OF    
  30.  * THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF THE AUTHOR HAS BEEN ADVISED   
  31.  * OF THE POSSIBILITY OF SUCH DAMAGE.                                         
  32.  *                                                                             
  33.  * THE AUTHOR SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT       
  34.  * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A    
  35.  * PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS"       
  36.  * BASIS, AND THE AUTHOR HAS NO OBLIGATIONS TO PROVIDE MAINTENANCE, SUPPORT,  
  37.  * UPDATES, ENHANCEMENTS, OR MODIFICATIONS.                                   
  38.  *                                                                             
  39.  *                                                                             
  40.  * [xsat for X-windows] Copyright 1992 by David A. Curry                                            
  41.  *                                                                             
  42.  * Permission to use, copy, modify, distribute, and sell this software and its 
  43.  * documentation for any purpose is hereby granted without fee, provided that  
  44.  * the above copyright notice appear in all copies and that both that copyright
  45.  * notice and this permission notice appear in supporting documentation.  The  
  46.  * author makes no representations about the suitability of this software for  
  47.  * any purpose.  It is provided "as is" without express or implied warranty.   
  48.  *                                                                             
  49.  * David A. Curry, N9MSW                                                       
  50.  * Purdue University                                                           
  51.  * Engineering Computer Network                                                
  52.  * 1285 Electrical Engineering Building                                        
  53.  * West Lafayette, IN 47907                                                    
  54.  * davy@ecn.purdue.edu                                                         
  55.  *                                                                             
  56.  * VersaTrack Copyright 1993, 1994 Siamack Navabpour. All Rights Reserved.
  57.  *
  58.  * Permission is hereby granted to copy, modify and distribute VersaTrack
  59.  * in whole, or in part, for educational, non-profit and non-commercial use
  60.  * only, free of charge or obligation, and without agreement, provided that
  61.  * all copyrights and restrictions noted herein are observed and followed, and
  62.  * additionally, that this and all other copyright notices listed herein
  63.  * appear unaltered in all copies and in all derived work.
  64.  *
  65.  * This notice shall not in any way void or supersede any of the other authors
  66.  * rights or privileges.
  67.  *
  68.  * VersaTrack IS PRESENTED FREE AND "AS IS", WITHOUT ANY WARRANTY OR SUPPORT.
  69.  * YOU USE IT AT YOUR OWN RISK. The author(s) shall not be liable for any
  70.  * direct, indirect, incidental, or consequential damage, loss of profits or
  71.  * other tangible or intangible losses or benefits, arising out of or related
  72.  * to its use. VersaTrack carries no warranty, explicit or implied, including
  73.  * but not limited to those of merchantablity and fitness for a particular
  74.  * purpose.
  75.  *
  76.  * Siamack Navabpour, 12342 Hunter's Chase Dr. Apt. 2114, Austin, TX 78729.
  77.  * sia@bga.com or sia@realtime.com.
  78.  */
  79.  
  80. /* THE FOLLOWING ADAPTED FROM Sattrack (the Unix version by Dr. M. Bester) */
  81.   
  82. #include <windows.h>
  83. #include <math.h>
  84.  
  85. #include "vstdefs.h"
  86. #include "vsttype.h"
  87. #include "constant.h"
  88. #include "vstextrn.h"
  89.  
  90. double
  91. absol( x, y, z)
  92. double x, y, z;
  93. {
  94.     return sqrt(x*x + y*y + z*z);
  95. }
  96.  
  97. double
  98. reduce(value, rangeMin, rangeMax)
  99. double value, rangeMin, rangeMax;
  100. {
  101.     double range, rangeFrac, fullRanges, retval;
  102.  
  103.     range     = rangeMax - rangeMin;
  104.     rangeFrac = (rangeMax - value) / range;
  105.  
  106.     modf(rangeFrac,&fullRanges);
  107.  
  108.     retval = value + fullRanges * range;
  109.  
  110.     if (retval > rangeMax)
  111.         retval -= range;
  112.  
  113.     return retval;
  114. }
  115.  
  116. static double
  117. getkep(ma, eccentricity, keperr)
  118. double ma, eccentricity;
  119. int *keperr;
  120. {
  121.     double eccAnomaly;                        /* eccentric anomaly */
  122.     double ta;
  123.     register double error;
  124.     int count;
  125.  
  126.     count = 5000;
  127.     eccAnomaly = ma;                          /* Initial guess     */
  128.     *keperr = 0;
  129.     do {
  130.         error = (eccAnomaly - eccentricity * sin(eccAnomaly) - ma) /
  131.                 (1.0 - eccentricity * cos(eccAnomaly));
  132.         eccAnomaly -= error;
  133.     } while (ABS(error) >= EPSILON && --count > 0);
  134.  
  135.     if (count <= 0)
  136.         *keperr = 1;
  137.  
  138.     if (ABS(eccAnomaly - PI) < EPSILON)
  139.         ta = PI;
  140.     else
  141.         ta = 2.0 * atan(sqrt((1.0 + eccentricity) /
  142.                             (1.0 - eccentricity)) * tan(eccAnomaly / 2.0));
  143.  
  144.     return  reduce (ta, (double)0.0, TWOPI);
  145. }
  146.  
  147. /*
  148.  * subsatpos: calculates sub-satellite point
  149.  *
  150.  *
  151.  *  EARTHRADIUS = equatorial radius of the Earth (m)
  152.  *                from "Astronomical Almanac", 1984, p. S7
  153.  *
  154.  *  F,FF,FFF    = flattening factors for the Earth (Hayford's spheroid)
  155.  *                from "Astronomical Almanac", 1984, p. S7
  156.  *  F1,2,3      = coefficients for the latitude correction formula
  157.  *                from "Explanatory Supplement to the Ephemeris", 1961, p.58
  158.  *
  159.  *  REFF0,1,2,3 = coefficient for the Earth's radius correction formula
  160.  *                from "Explanatory Supplement to the Ephemeris", 1961, p.58
  161.  *
  162.  *  effEarthRad = effective radius of the Earth at the location of the
  163.  *                sub-satellite point with respect to sea level (m)
  164.  *
  165.  *  for reference see:  a) "Methods of Experimental Physics", Vol. 12C, p.308
  166.  *                      b) "Explanatory Supplement to the Ephemeris",1961, p.58
  167.  */
  168.  
  169. static void
  170. ssp(tp, satX, satY, satZ, time, latitude, longitude, height, avis)
  171. track_t *tp;
  172. double satX, satY, satZ, time;
  173. float *latitude, *longitude;
  174. double *height, *avis;
  175. {
  176.     double x, y, z, r, b2, b4, b6, lgt, effrad;
  177.  
  178.     r   = absol(satX, satY, satZ);
  179.  
  180.     x   = satX / r;
  181.     y   = satY / r;
  182.     z   = satZ / r;
  183.  
  184.     if (z >  1.0) z =  1.0;
  185.     if (z < -1.0) z = -1.0;
  186.  
  187.     *latitude   = (float) asin(z);
  188.  
  189.     b2          = 2.0 * (*latitude);
  190.     b4          = b2 + b2;
  191.     b6          = b2 + b4;
  192.  
  193.     *latitude  += (float) ((F1 * sin(b2)) + (F2 * sin(b4)) + (F3 * sin(b6)));
  194.  
  195. #ifdef SIMPLER_MODEL
  196.     lgt         = TWOPI * ((time - tp->sidday) * SIDSOLAR + tp->sidtime)
  197.                   - atan2(y,x);
  198. #else
  199.     lgt         = tp->gastime - atan2(y,x);
  200. #endif
  201.  
  202.     *longitude  = (float) reduce(lgt, -PI, PI);
  203.  
  204.     effrad      = (REFF0 + REFF1 * cos(b2) + REFF2 * cos(b4) + REFF3 * cos(b6))
  205.                   * EARTHRADIUS;
  206.  
  207.     *height     = r - effrad;
  208.     
  209.     /* XXX FIX - sia: set crash flag if height <= 0.0 */
  210.     
  211.     *avis       = acos(effrad/r);   /* geocenteric angle of visibility */
  212. }
  213.  
  214. /*
  215.  * precess: calculates precession of the satellite's orbit
  216.  */
  217.  
  218. static void
  219. precess(semiMajorAxis, eccentricity, inclination, RAANPrec, perigeePrec)
  220. double semiMajorAxis, eccentricity, inclination;
  221. double *RAANPrec, *perigeePrec;
  222. {
  223.     double factor;
  224.  
  225.     /*
  226.      * See the Satellite Experimenter's Handbook, pp. 11-11,11-12,11-13.
  227.      * The 3rd mulitiplicative term in the formula is reduced to zero
  228.      * when orbit is near circluar. i.e., eccentricity almost = 0.
  229.      * Equation 11.16.
  230.      */
  231.      
  232.     factor       = pow(EARTHRADIUS / semiMajorAxis, 3.5)
  233.                    / SQR(1.0 - SQR(eccentricity)) * CDR;
  234.  
  235.     *RAANPrec    = 9.95 * factor * cos(inclination);
  236.  
  237.     /*
  238.      * See Satellite Experimenter's Handbook, p. 11-9. Equation 11.13a
  239.      */
  240.     *perigeePrec = 4.97 * factor * ( 5.0 * SQR( cos(inclination) ) - 1.0 );
  241. }
  242.  
  243. /*
  244.  * satpos : calculates satellite position and velocity in the mean
  245.  *          equatorial coordinate system (RA, Dec)
  246.  */
  247.  
  248. static int
  249. satpos(epochTime, epochRAAN, epochArgPerigee, semiMajorAxis,
  250.         inclination, eccentricity, RAANPrec, perigeePrec, time, trueAnomaly,
  251.         curRaan, argPerigee, X, Y, Z, radius, vX, vY, vZ)
  252.  
  253. double epochTime, epochRAAN, epochArgPerigee;
  254. double semiMajorAxis, inclination, eccentricity;
  255. double RAANPrec, perigeePrec, time, trueAnomaly;
  256. double *curRaan, *argPerigee, *X, *Y, *Z, *radius, *vX, *vY, *vZ;
  257. {
  258.     double RAAN;
  259.     double Xw, Yw, vXw, vYw;                    /* in orbital plane */
  260.     double tmp;
  261.     double Px, Py, Pz, Qx, Qy, Qz;              /* Escobal transformation #31 */
  262.     double cosArgPerigee, sinArgPerigee;
  263.     double cosRAAN, sinRAAN, cosInclination, sinInclination;
  264.  
  265.     /* See Satellite Experimenter's Handbook p. 11-6. */
  266.     
  267.     *radius = semiMajorAxis * (1.0 - SQR(eccentricity))
  268.               / (1.0 + (eccentricity * cos(trueAnomaly)));
  269.  
  270.     if (*radius < EARTHRADIUS) {
  271. #ifdef _DEBUG_
  272.         diag("satellite has crashed!\n");
  273. #endif /* _DEBUG_ */
  274.         usermsg("Satellite not in orbit! Elementset may be old or incorrect");
  275.         return 0;
  276.     }
  277.  
  278.     Xw              = *radius * cos(trueAnomaly);
  279.     Yw              = *radius * sin(trueAnomaly);
  280.     tmp             = sqrt(GMSGP / (semiMajorAxis * (1.0 - SQR(eccentricity))));
  281.     vXw             = -tmp * sin(trueAnomaly);
  282.     vYw             = tmp * (cos(trueAnomaly) + eccentricity);
  283.  
  284.     *argPerigee     = epochArgPerigee + (time - epochTime) * perigeePrec;
  285.     *curRaan = RAAN = epochRAAN - (time - epochTime) * RAANPrec;
  286.     cosRAAN         = cos(RAAN);
  287.     sinRAAN         = sin(RAAN);
  288.     cosArgPerigee   = cos(*argPerigee);
  289.     sinArgPerigee   = sin(*argPerigee);
  290.     cosInclination  = cos(inclination);
  291.     sinInclination  = sin(inclination);
  292.     
  293.     Px  =  cosArgPerigee * cosRAAN - sinArgPerigee * sinRAAN * cosInclination;
  294.     Py  =  cosArgPerigee * sinRAAN + sinArgPerigee * cosRAAN * cosInclination;
  295.     Pz  =  sinArgPerigee * sinInclination;
  296.  
  297.     Qx  = -sinArgPerigee * cosRAAN - cosArgPerigee * sinRAAN * cosInclination;
  298.     Qy  = -sinArgPerigee * sinRAAN + cosArgPerigee * cosRAAN * cosInclination;
  299.     Qz  =  cosArgPerigee * sinInclination;
  300.  
  301.     *X  =  Px*Xw + Qx*Yw;                       /* Escobal transformation #31 */
  302.     *Y  =  Py*Xw + Qy*Yw;
  303.     *Z  =  Pz*Xw + Qz*Yw;
  304.  
  305.     *vX =  Px*vXw + Qx*vYw;                     /* satellite velocity         */
  306.     *vY =  Py*vXw + Qy*vYw;
  307.     *vZ =  Pz*vXw + Qz*vYw;
  308.     return 1;
  309. }
  310.  
  311. /*
  312.  * sitepos: computes the site postion and velocity in the mean
  313.  *          equatorial coordinate system. The matrix 'siteMat' is
  314.  *          used by geototopo to convert geocentric coordinates
  315.  *          to topocentric (observer-centered) coordinates.
  316.  */
  317.  
  318. static void
  319. sitepos(tp, siteLat, siteLong, siteAlt, siteTime,
  320.                 siteX, siteY, siteZ, siteVX, siteVY)
  321. track_t *tp;
  322. double siteLat, siteLong, siteAlt, siteTime;
  323. double *siteX, *siteY, *siteZ, *siteVX, *siteVY;
  324. {
  325.     double sb2, sb4, sb6, siteRA, cosRA, sinRA, effSiteRad;
  326.     double G1, G2, cosLat, sinLat;
  327.  
  328.     sb2        = 2.0 * siteLat;
  329.     sb4        = sb2 + sb2;
  330.     sb6        = sb2 + sb4;
  331.  
  332. #ifdef SIMPLER_MODEL
  333.     siteLat   += F1 * sin(sb2) + F2 * sin(sb4) + F3 * sin(sb6);
  334. #endif    
  335.     cosLat     = cos(siteLat);
  336.     sinLat     = sin(siteLat);
  337.  
  338.     effSiteRad = (REFF0 + REFF1 * cos(sb2) + REFF2 * cos(sb4) + 
  339.                           REFF3 * cos(sb6)) * EARTHRADIUS + siteAlt;
  340.  
  341.     G1         = effSiteRad * cosLat;
  342.     G2         = effSiteRad * sinLat;
  343. #ifdef SIMPLER_MODEL
  344.     siteRA  =  TWOPI * ((siteTime - tp->sidday) * SIDSOLAR + tp->sidtime) - siteLong;
  345. #else
  346.     siteRA  =  tp->lastime;
  347. #endif
  348.     cosRA   =  cos(siteRA);
  349.     sinRA   =  sin(siteRA);
  350.  
  351.     *siteX  =  G1 * cosRA;
  352.     *siteY  =  G1 * sinRA;
  353.     *siteZ  =  G2;
  354.  
  355.     *siteVX = -SIDRATE * (*siteY);
  356.     *siteVY =  SIDRATE * (*siteX);
  357.  
  358.     tp->sitemat[0][0] =  sinLat * cosRA;
  359.     tp->sitemat[0][1] =  sinLat * sinRA;
  360.     tp->sitemat[0][2] = -cosLat;
  361.     tp->sitemat[1][0] = -sinRA;
  362.     tp->sitemat[1][1] =  cosRA;
  363.     tp->sitemat[1][2] =  0.0;
  364.     tp->sitemat[2][0] =  cosRA * cosLat;
  365.     tp->sitemat[2][1] =  sinRA * cosLat;
  366.     tp->sitemat[2][2] =  sinLat;
  367. }
  368.  
  369. /*
  370.  * satange: calculates satellite range
  371.  */
  372.  
  373. #ifdef RANGE_USEROUTINE
  374.  
  375. void satrange(siteX,siteY,siteZ,siteVX,siteVY,
  376.               satX,satY,satZ,satVX,satVY,satVZ,range,rangeRate)
  377. double siteX, siteY, siteZ, siteVX, siteVY;
  378. double satX, satY, satZ, satVX, satVY, satVZ;
  379. double *range, *rangeRate;
  380. {
  381.     double dX, dY, dZ;
  382.  
  383.     dX = satX - siteX;
  384.     dY = satY - siteY;
  385.     dZ = satZ - siteZ;
  386.  
  387.     *range     = absol(dX ,dY dZ);    
  388.     *rangeRate = ((satVX-siteVX)*dX + (satVY-siteVY)*dY + satVZ*dZ) / *range;
  389. }
  390. #endif /* RANGE_USEROUTINE */
  391.  
  392. /*
  393.  * geototopo: converts from geocentric mean equatorial coordinates to
  394.  *            topocentric coordinates
  395.  */
  396.  
  397. static void
  398. geototopo(satX, satY, satZ, siteX, siteY, siteZ, siteMatT, X, Y, Z)
  399. double satX, satY, satZ, siteX, siteY, siteZ;
  400. double *X, *Y, *Z;
  401. double siteMatT[3][3];
  402. {
  403.     satX -= siteX;
  404.     satY -= siteY;
  405.     satZ -= siteZ;
  406.     *X    = siteMatT[0][0]* satX + siteMatT[0][1]* satY + siteMatT[0][2]* satZ;
  407.     *Y    = siteMatT[1][0]* satX + siteMatT[1][1]* satY + siteMatT[1][2]* satZ;
  408.     *Z    = siteMatT[2][0]* satX + siteMatT[2][1]* satY + siteMatT[2][2]* satZ;
  409. }
  410.  
  411. /*
  412.  * ae: calculates azimuth and elevation
  413.  */
  414. static void
  415. ae(satX, satY, satZ, siteX, siteY, siteZ, siteMatA, azimuth, elevation)
  416. double satX, satY, satZ, siteX, siteY, siteZ;
  417. double siteMatA[3][3];
  418. float *azimuth, *elevation;
  419. {
  420.     double  x, y, z, r;
  421.  
  422.     geototopo(satX, satY, satZ, siteX, siteY, siteZ, siteMatA, &x, &y, &z);
  423.     
  424.     r  = absol(x ,y ,z);
  425.     
  426.     x /= r;
  427.     y /= r;
  428.     z /= r;
  429.  
  430.     if (x == ZERO)
  431.         *azimuth = (y > ZERO) ? (float) HALFPI : (float)  THREEHALFPI;
  432.     else        
  433.         *azimuth = (float) (PI - atan2(y, x));
  434.  
  435.     if (z >  1.0) z =  1.0;      /* protect against out-of-range problems */
  436.     if (z < -1.0) z = -1.0;
  437.  
  438.     *elevation = (float) asin(z);
  439. }
  440.  
  441. /*
  442.  * isinsun: calculates the eclipses of the satellite
  443.  *
  444.  * this function assumes that the Sun is point-like rather than extended
  445.  * and that the Earth is perfectly spherical
  446.  * XXX sia - DOUBLE CHECK FOLLOWING.
  447.  */
  448.  
  449. static int
  450. insun(tp, satX, satY, satZ, siteX, siteY, siteZ, eclTime, satElevation)
  451. track_t *tp;
  452. double satX, satY, satZ, siteX, siteY, siteZ, eclTime;
  453. float satElevation;
  454. {
  455.     double meanAnomaly, trueAnomaly;
  456.     double sunX, sunY, sunZ, sunVX, sunVY, sunVZ, sunRadius, sunDist;
  457.     double satPara, satPerp, satDist, alpha, beta, fdummy;
  458.     int err;
  459.  
  460.     meanAnomaly  = tp->sunmeananomaly +  ((eclTime - tp->sunepochtime) *
  461.                    tp->sunmeanmotion * TWOPI);
  462.     trueAnomaly  = getkep(meanAnomaly, tp->suneccentricity, &err);
  463.     
  464.     satpos(tp->sunepochtime, tp->sunRAAN, tp->sunargperigee,
  465.         SUNSEMIMAJORAXIS, tp->suninclination, tp->suneccentricity, 0.0, 0.0,
  466.         eclTime, trueAnomaly,  &fdummy, &fdummy, &sunX, &sunY, &sunZ, &sunRadius,
  467.         &sunVX, &sunVY, &sunVZ);
  468.  
  469.     ae(sunX, sunY, sunZ, siteX, siteY, siteZ, tp->sitemat,
  470.         &tp->sunazimuth, &tp->sunelevation);
  471.  
  472.     sunDist = absol(sunX, sunY, sunZ);
  473.     satDist = absol(satX, satY, satZ);
  474.  
  475.     satPara = (satX*sunX + satY*sunY + satZ*sunZ) / sunDist;
  476.     satPerp = sqrt(SQR(satDist) - SQR(satPara));
  477.  
  478.     alpha = asin(EARTHRADIUS / sunDist);
  479.     beta  = atan(satPerp / (satPara + sunDist));
  480.  
  481.     /* case 1: satellite is on side of the Earth facing the Sun */
  482.  
  483.     if (satPara > 0.0) {
  484. #ifdef ORIGINAL_CODE
  485.         if (satElevation < 0.0)
  486.             return(DAY);
  487.   
  488.         if (satElevation >= 0.0 && tp->sunelevation <= TWILIGHT)
  489.             return(VISIBLE);
  490.         else
  491. #endif
  492.             return (int) 'Y'; /* in sunlight */
  493.     }
  494.  
  495.     /* case 2: satellite is on side of the Earth facing away from the Sun */
  496.  
  497.     if (beta >= alpha) {   /* satellite can see the Sun */
  498.         return (int) 'Y';
  499. #ifdef ORIGINAL_CODE
  500.         if (satElevation >= 0.0 && tp->sunelevation <= TWILIGHT)
  501.             return(VISIBLE);
  502.         else
  503.             return(DAY);
  504. #endif
  505.     }
  506.     else                  /* satellite is in the Earth's shadow */
  507.         return (int) 'N'; /* not in sun */
  508. }
  509.  
  510. /*
  511.  * confss: initializes the Sun's Keplerian elements for a given
  512.  *         epoch. Formulas are from "Explanatory Supplement to the
  513.  *         Ephemeris, HMSO, 1974".
  514.  *         It also initializes the sidereal reference when the simpler
  515.  *         model is used.
  516.  */
  517.  
  518. void confss(tp, initTime)
  519. track_t *tp;
  520. double initTime;
  521. {
  522.     double T, T2, T3, omega;
  523.     int n;
  524.  
  525.     if (fabs(initTime - tp->lastsuntime) <= SUNUPDATEINT)
  526.         return;
  527.  
  528.     tp->lastsuntime = initTime;
  529.  
  530.     /* time argument in Julian centuries since 1900.0 */
  531.     /* see Astronomical Almanac, 1983, page B6 */
  532.     
  533.     /* for nutation correction see Explanatory Supplement, 1974, page 44 */
  534.  
  535.     tp->sunRAAN     = 0.0;
  536.     tp->sunepochtime= initTime;
  537.     tp->sidday      = floor(initTime);
  538.  
  539.     T               = (initTime - 0.5) / 36525.0;
  540.     T2              = T*T;
  541.     T3              = T2*T;
  542.  
  543.     tp->suneccentricity = 0.01675104 - 4.180e-5*T - 1.260e-7*T2;
  544.     tp->sunargperigee   = (281.220844 + 1.719175*T + 4.528e-4*T2 +
  545.                            3.3333e-6*T3) * CDR;
  546.  
  547.     tp->sunmeanmotion   = 1.0 / (365.24219878 - 6.134e-6*T);
  548.  
  549.     omega           = (259.183275 - 1934.142008*T + 2.078e-3*T2 +
  550.                        2.22e-6*T3) * CDR;
  551.     n               = (int) (omega / TWOPI);
  552.     omega          -= (double) n * TWOPI;
  553.  
  554.     tp->suninclination  = ((23.452294 - 0.0130125*T - 1.64e-6*T2 + 5.03e-7*T3) +
  555.                           2.558333e-3*cos(omega)) * CDR;       /* page 41 */
  556. #ifdef SIMPLER_MODEL
  557.     tp->sidtime     = (6.646065556 + 2400.051262*T + 2.580556e-5*T2) / 24.0;
  558.     tp->sidtime    -= floor(tp->sidtime);
  559. #else
  560.     /* sia - use local apparent sidreal computed in getJulianDate */
  561. #endif
  562.  
  563.     /* for ephemeris of the Sun see Explanatory Supplement, 1974, page 98 */
  564.  
  565.     tp->sunmeananomaly  = (358.475833 + 35999.04975*T - 1.5e-4*T2 -
  566.                            3.3333e-6*T3) * CDR;
  567.     n                   = (int) (tp->sunmeananomaly / TWOPI);
  568.     tp->sunmeananomaly -= (double) n * TWOPI;
  569.     tp->suntrueanomaly  = getkep(tp->sunmeananomaly, tp->suneccentricity, &n);
  570.  
  571.     tp->sundistance     = SUNSEMIMAJORAXIS * (1.0 - SQR(tp->suneccentricity))
  572.                       / (1.0 + tp->suneccentricity * cos(tp->suntrueanomaly));
  573. }
  574.  
  575. /*
  576.  * Theoretical free space signal attenuation. The result is *NOT* a good
  577.  * indication of the link budget. It does not take into consideration
  578.  * the effects of atmospheric composition (signal scattering and absorbtion),
  579.  * ionspheric conditions (loss due to changes in polarization, and effects
  580.  * of high energy particles, etc.), antenna apprature and gain parameters,
  581.  * etc. etc. SN101094.
  582.  */
  583. static float
  584. attenuation(range, freq)
  585. double range, freq;
  586. {
  587.     return  (float) (20.0 * log10 ( freq * ((FOURPI * range) / CVAC) ));
  588. }
  589.  
  590. static int
  591. ta(timeArg, tp)     /* True Anomaly */
  592. double timeArg;
  593. track_t *tp;
  594. {
  595.     double deltaT, fdummy, curMotion, curOrbit;
  596.     int error;
  597.  
  598.     deltaT            = timeArg - tp->epochday;
  599.     curMotion         = tp->epochmeanmotion + deltaT * tp->decayrate;
  600.     tp->semimajoraxis = KEPLER * exp(2.0 * log(MPD / curMotion) / 3.0);
  601.     curOrbit          = tp->reforbit + deltaT * curMotion;
  602.     tp->orbitnum      = (long) curOrbit;
  603.     tp->meananomaly   = (curOrbit - (double) tp->orbitnum) * TWOPI;
  604.     tp->trueanomaly   = getkep(tp->meananomaly, tp->eccentricity, &error);
  605.     tp->orbitfrac     = (int) (modf(curOrbit, &fdummy) * 100.0);
  606.     tp->orbitnum     += 1;
  607.     return error;
  608. }
  609.  
  610.  
  611. void satphase(tp)
  612. track_t *tp;
  613. {
  614.     int n, phase;
  615.     double dphase;
  616.  
  617.     dphase = (tp->meananomaly / TWOPI * tp->maxphase + tp->perigeephase);
  618.     phase  = (int) dphase;
  619.  
  620.     n = 0;
  621.     while (phase < 0 && ++n < 1000)
  622.         phase += (int) tp->maxphase;
  623.  
  624.     n = 0;
  625.     while (phase >= (int) tp->maxphase && ++n < 1000)
  626.         phase -=  (int) tp->maxphase;
  627.  
  628.     tp->phase = phase;
  629. }
  630.  
  631. /*
  632.  * satmode: gets satellite mode from phase
  633.  */
  634.  
  635. void satmode(sp, phase, modestrp)
  636. satellite_t *sp;
  637. int phase;
  638. char *modestrp;
  639. {
  640.     mode_t *mp;
  641.     int curMode;
  642.  
  643.     if (sp && (mp = sp->s_modep) ) {
  644.         for (curMode = 0; curMode < mp->nmodes; curMode++) {
  645.             if ( ( phase >= mp->mds[curMode].minPhase  &&
  646.                 phase < mp->mds[curMode].maxPhase)  ||
  647.                 ( (mp->mds[curMode].minPhase > mp->mds[curMode].maxPhase) &&
  648.                 ( phase >= mp->mds[curMode].minPhase ||
  649.                 phase < mp->mds[curMode].maxPhase) ) ) {
  650.  
  651.                 strncpy(modestrp, mp->mds[curMode].modeStr,2);
  652.                 return;
  653.             }
  654.         }
  655.     }
  656.     strcpy(modestrp,"N/A");
  657. }
  658.  
  659.  
  660. void
  661. predict_init(flag, sp, ptime)
  662. int flag;
  663. select_t *sp;
  664. double ptime;
  665. {
  666.     track_t *tp;
  667.     result_t *rp;
  668.  
  669.     tp = sp->sl_tp;
  670.     rp = sp->sl_rp;
  671.  
  672.     params_init(tp->satp, tp);
  673.     confss(tp, ptime);    
  674.     siderealtime(tp, ptime);
  675.     
  676.     tp->semimajoraxis = KEPLER * exp(2.0 * log(MPD / tp->epochmeanmotion) / 3.0);
  677.     tp->reforbit = (double) tp->epochorbitnum + tp->epochmeananomaly / TWOPI;
  678.   
  679.     predict(flag, sp, ptime);
  680. }
  681.  
  682.  
  683. int
  684. predict(flag, sp, ptime)
  685. int flag;
  686. select_t *sp;
  687. double ptime;
  688. {
  689.     site_t *sitep;
  690.     satellite_t *satp;
  691.     track_t *tp;
  692.     result_t *rp;
  693.     double siteX, siteY, siteZ, siteVX, siteVY, argperigee, curRaan;
  694.     double satX, satY, satZ, satVX, satVY, satVZ, satRadius, dX,dY,dZ;
  695.     double avis;
  696.     extern void satposSGP4(select_t *, double, double *, double *, double *,
  697.         double *,  double *, double *, double *);
  698.     extern void SquintAngle(select_t *, double, double, double, double);
  699.  
  700.     tp = sp->sl_tp;
  701.     rp = sp->sl_rp;
  702.     sitep = tp->sitep;
  703.     satp = tp->satp;
  704.     
  705.     confss(tp, ptime);
  706.  
  707.     siderealtime(tp, ptime);
  708.  
  709.     if (fabs(tp->lastprectime - ptime) > PRECUPDATEINT) {
  710.         precess(tp->semimajoraxis, tp->eccentricity, tp->inclination,
  711.                        &tp->RAANprec, &tp->perigeeprec);
  712.         tp->lastprectime = ptime;
  713.     }
  714.  
  715.     sitepos(tp, sitep->c_lat, sitep->c_lng, sitep->c_alt, ptime,
  716.         &siteX, &siteY, &siteZ, &siteVX, &siteVY);
  717.  
  718. #ifdef RARELYIFEVERATALL
  719.     if (ta(ptime, tp))
  720.         return 0;
  721.  
  722.     if (!satpos(tp->epochday, tp->epochRAAN, tp->epochargperigee,
  723.         tp->semimajoraxis,  tp->inclination, tp->eccentricity,
  724.         tp->RAANprec, tp->perigeeprec, ptime, tp->trueanomaly, &curRaan,
  725.         &argperigee, &satX, &satY, &satZ, &satRadius, &satVX, &satVY, &satVZ))
  726.                 return 0;
  727.  
  728.     ae(satX, satY, satZ, siteX, siteY, siteZ, tp->sitemat, &rp->r_azimuth,
  729.         &rp->r_elevation);
  730.  
  731.     printf("MEAN  X:  %.3lf   Y: %.3lf  Z: %lf\n", satX, satY, satZ);
  732.     printf("      AZ: %.3lf  EL: %.3lf\n", rp->r_azimuth, rp->r_elevation);
  733.  
  734.     satposSGP4(sp, ptime, &satX, &satY, &satZ, &satRadius, &satVX, &satVY, &satVZ);
  735.  
  736.     ae(satX, satY, satZ, siteX, siteY, siteZ, tp->sitemat, &rp->r_azimuth,
  737.         &rp->r_elevation);
  738.  
  739.     printf("SGP4  X:  %.3lf   Y: %.3lf  Z: %lf\n", satX, satY, satZ);
  740.     printf("      AZ: %.3lf  EL: %.3lf\n", rp->r_azimuth, rp->r_elevation);
  741.  
  742. #else
  743.     if (flag == TLEMEAN) {
  744.         if (ta(ptime, tp))
  745.             return 0;
  746.  
  747.         if (!satpos(tp->epochday, tp->epochRAAN, tp->epochargperigee,
  748.             tp->semimajoraxis,  tp->inclination, tp->eccentricity,
  749.             tp->RAANprec, tp->perigeeprec, ptime, tp->trueanomaly,
  750.             &curRaan, &argperigee, &satX, &satY, &satZ, &satRadius,
  751.             &satVX, &satVY, &satVZ))
  752.                 return 0;
  753.             SquintAngle(sp, tp->trueanomaly, curRaan, argperigee, satRadius);
  754.     }
  755.     else if (flag == NORADSGP4)
  756.         satposSGP4(sp, ptime, &satX, &satY, &satZ, &satRadius, &satVX, &satVY, &satVZ);
  757.     else {
  758.         ; /* deep-space (!) prop model. */
  759.     }
  760.     ae(satX, satY, satZ, siteX, siteY, siteZ, tp->sitemat, &rp->r_azimuth,
  761.         &rp->r_elevation);
  762.  
  763. #endif /* _DEBUG_ */
  764.  
  765. #ifdef RANGE_USEROUTINE
  766.     satrange(siteX, siteY, siteZ, siteVX, siteVY, satX, satY, satZ,
  767.         satVX, satVY, satVZ, &rp->r_range, &rp->r_rangerate);
  768. #else   
  769.     dX = satX - siteX;
  770.     dY = satY - siteY;
  771.     dZ = satZ - siteZ;
  772.  
  773.     rp->r_range    = absol(dX, dY, dZ);    
  774.     rp->r_rangerate = ((satVX-siteVX)*dX + (satVY-siteVY)*dY + satVZ*dZ) /
  775.             rp->r_range;
  776. #endif /* RANGE_USEROUTINE */
  777.  
  778.     ssp(tp, satX, satY, satZ, ptime, &rp->r_lat, &rp->r_long,
  779.         &rp->r_height, &avis);
  780.     
  781.     rp->r_doppler = (float) ((- tp->beacon * rp->r_rangerate) / CVAC); /* in HZ */
  782.  
  783.     rp->r_pathloss = attenuation(rp->r_range, tp->beacon);
  784.  
  785.     satphase(tp);
  786.  
  787.     satmode(satp, rp->r_phase = tp->phase, rp->r_modestr); 
  788.  
  789.     rp->r_insun = (char) insun(tp, satX, satY, satZ, siteX, siteY, siteZ,
  790.         ptime, rp->r_elevation);
  791.  
  792.     rp->r_av = avis * CRD;
  793.     rp->r_orbitnum = tp->orbitnum;
  794.     rp->r_lat *= (float) (CRD);
  795.     rp->r_long *= (float) (CRD);
  796.     rp->r_pelev = rp->r_elevation;
  797.     rp->r_elevation *= (float) (CRD);
  798.     rp->r_azimuth *= (float)(CRD);
  799.     rp->r_time = ptime;
  800.     rp->r_rising = (rp->r_elevation >= rp->r_pelev) ? TRUE : FALSE;
  801.     return 1;
  802. }
  803.  
  804.  
  805. static void
  806. params_init(satp, trackp)
  807. satellite_t *satp;
  808. track_t *trackp;
  809. {        
  810.     mode_t *mp;
  811.  
  812.     trackp->satname          = satp->s_name;
  813.     trackp->satnum           = satp->s_satnum;
  814.     trackp->elementset       = satp->s_elementset;
  815.     trackp->epochday         = satp->s_epochtime;
  816.     trackp->inclination      = satp->s_inclination * CDR;
  817.     trackp->epochRAAN        = satp->s_raan * CDR;
  818.     trackp->eccentricity     = satp->s_eccentricity;
  819.     trackp->epochargperigee  = satp->s_argperigee * CDR;
  820.     trackp->epochmeananomaly = satp->s_meananomaly * CDR;
  821.     trackp->epochmeanmotion  = satp->s_meanmotion;
  822.     trackp->decayrate        = satp->s_decayrate;
  823.     trackp->epochorbitnum    = satp->s_orbitnum;
  824.     trackp->flags            = 0;
  825.     trackp->lastjdate        = ZERO;
  826.     trackp->lastsuntime      = ZERO;
  827.     trackp->lastprectime     = ZERO;
  828.  
  829.     if ((fabs(satp->s_meanmotion - 1.0) < 0.05) &&
  830.         (fabs(satp->s_inclination) < 1.0))          /* XXX FIX - sia */
  831.         satp->s_flags |= SF_GEOSTAT;
  832.  
  833.     if (satp && (mp=satp->s_modep)) {
  834.         trackp->beacon       = mp->beacon;
  835.         trackp->perigeephase = mp->perigeePhase;
  836.         trackp->maxphase     = mp->maxPhase;
  837.     }
  838.     else {
  839.         trackp->beacon       = BEACON * 1.0e6;
  840.         trackp->perigeephase = 0.9;
  841.         trackp->maxphase     = MAXPHASE;
  842.     }
  843. }
  844.  
  845. void
  846. printMode(fp, sp, phase)
  847. FILE *fp;
  848. satellite_t *sp;
  849. int phase;
  850. {
  851.     int curMode;
  852.     mode_t *mp;
  853.  
  854.     if (sp && (mp=sp->s_modep)) {     
  855.         for (curMode = 0; curMode < mp->nmodes; curMode++) {
  856.             if ((phase >= mp->mds[curMode].minPhase
  857.                 && phase < mp->mds[curMode].maxPhase)
  858.                 || ((mp->mds[curMode].minPhase > mp->mds[curMode].maxPhase)
  859.                 && (phase >= mp->mds[curMode].minPhase
  860.                 || phase < mp->mds[curMode].maxPhase))) {
  861.                 
  862.                 fprintf(fp," %2s",mp->mds[curMode].modeStr);
  863.             }
  864.         }
  865.     }
  866.     else
  867.         fprintf(fp,"  ");
  868. }
  869.  
  870.  
  871.